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Fourier methods are fundamental tools to analyze random fields. Statistical structures of homoge- 
neous Gaussian random fields are completely characterized by the power spectrum. In non-Gaussian 
random fields, polyspectra, higher-order counterparts of the power spectrum, are usually considered 
to characterize statistical information contained in the fields. However, it is not trivial how the 
Fourier modes are distributed in general non-Gaussian fields. In this paper, distribution functions 
of Fourier modes are directly considered and their explicit relations to the polyspectra are given. Un- 
der the condition that any of the polyspectra does not diverge, the distribution function is expanded 
by dimensionless combinations of polyspectra and a total volume in which the Fourier transforms are 
performed. The expression up to second order is generally given, and higher-order results are also 
derived in some cases. A concept of iV-point distribution function of Fourier modes are introduced 
and explicitly calculated. Among them, the one-point distribution function is completely given in 
a closed form up to arbitrary order. As an application, statistics of Fourier phases are explored in 
detail. A number of aspects regarding statistical properties of phases are found. It is clarified, for 
the first time, how phase correlations arise in general non-Gaussian fields. Some of our analytic 
results are tested against numerical realizations of non-Gaussian fields, showing good agreements. 

PACS numbers: 02.30.Nw, 05.40.-a, 98.65.Dx 



I. INTRODUCTION 



The Fourier analysis has long been one of the most important methods in various fields of study. Among them, 
randomly varying data are efficiently analyzed by Fourier transform methods, or spectral methods. In random 
processes as functions of time, the Fourier transform plays one of the central roles in unveiling the statistical nature 
of the process 0, 0> E ■ A random process is represented by a random function defined over 1-dimensional space. 
Generalizations to multi-dimensional spaces are straightforward, and random functions defined over some space are 
called random fields 0, Q . A random field is called homogeneous when all the statistical properties are translationally 
invariant and independent on a position in space. In a one-dimensional case, such as a random process, homogeneous 
random fields are also called stationary. Fourier analyses of random processes or random fields are especially useful 
when they are stationary or homogeneous. 

There is a class of random fields which are called Gaussian random fields. When a field arises from the superposition 
of a large number of independent random effects, the resulting field is a Gaussian random field under very weak 
conditions, by virtue of the so-called central limit theorem. There are many situations encountered in physical 
and engineering problems that random fields are accurately or approximately Gaussian. Statistical properties of a 
homogeneous Gaussian random field is completely specified by a correlation function, or equivalently, its Fourier 
transform, a power spectrum. Therefore, Gaussian random fields are analytically easier to treat than other random 
fields, and many aspects of Gaussian random fields have been investigated before |5j. 

In cosmological physics, Gaussian random fields are very important, since the initial density fluctuations are thought 
to be, at least approximately, a Gaussian random field. In inflationary models, the curvature perturbations generated 
by quantum fluctuations in the very early universe would yield a Gaussian random density field. Statistics of Gaussian 
random fields in cosmological contexts have been explored in detail 0, • 

However, non-Gaussianity naturally appears in reality, especially when nonlinear dynamics are involved. While lin- 
ear evolution of an initial density field in cosmology keeps the Gaussianity, nonlinear evolution raises non-Gaussianity 
in the density field 9]. Moreover, whether or not the initial density field is a Gaussian random field is an important 
question to investigate the beginning of the universe, and understanding non-Gaussian random fields has a great 
importance. 

The non-Gaussian features in Fourier analyses are usually characterized by polyspectra, i.e., higher-order counter- 
parts of the power spectra. For a given set of random variables, a set of all moments of the variables carries complete 
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information of the statistical distribution. However, it is still not trivial to understand the statistical properties 
of the Fourier modes, even when a hierarchy of all the polyspectra are known. For example, statistical behaviors of 
Fourier phases in non-Gaussian random fields has long been studied empirically, lackin g any w ell-established statistical 
framework to understand information contained in the Fourier phases [l(J, [HI U4> UJa LL-J uM ■ 

The statistical structures of Fourier phases are recently unveiled by the present author, properly calculating the joint 
distribution function of phases 0|. The analytical results are tested against numerical simulations of cosmological 
structure formation 01 > an d applied to an analysis of observational data • The analytical calculations involve an 
expansion by a certain parameter, and only lowest-order terms were discussed in the previous work. 

In this paper, the previous analytical calculations are extended to give a framework for studying statistics of Fourier 
modes in general. The statistical structures of not only Fourier phases but also Fourier moduli are simultaneously 
considered. We give a general joint distribution function of Fourier modes in general non-Gaussian fields. The 
joint distribution function of Fourier modes has complete information of statistical properties. Since a set of all the 
polyspectra also has the complete information, the joint distribution function can, in principle, be expressible by a 
set of all the polyspectra. Such relation has not been known, and it is surely difficult to obtain the complete relation. 
However, the distribution function is found to inevitably depend on a total volume in which Fourier transforms are 
applied. Some dimensionless combinations of a polyspectrum and a total volume, each of which is a dimensional 
quantity, can actually be arbitrary small for a large volume limit. Therefore, it is natural to expand the distribution 
function by some powers of a volume. The expression is found to be expanded by an inverse of square root of the 
volume as a small parameter. The way to calculate the distribution function up to arbitrary order is established in this 
paper. The general expression up to second order is actually given. More general properties involving higher-order 
effects are also obtained and discussed in some cases. Remarkably, an analytic expression of one-point distribution 
function of a Fourier mode is obtained to arbitrary order in a closed form. 

This paper is organized as follows. In Sec. [OJ our notations of Fourier transforms are defined, and the distribution 
function of Fourier modes in Gaussian random fields are reviewed. In Sec. lIIII methods to derive distribution function 
in general non-Gaussian fields are detailed, and an explicit expression up to second order is obtained. In Sec. IIVI TV- 
point distribution functions of Fourier modes are introduced, and calculated in several cases. An analytic expression 
of the one-point distribution function is completely given in a closed form. General expression of the two-, three-, 
four-point distribution functions are given up to second order. In Sec. |Vj statistical structures of Fourier phases are 
revealed. ./V-point distribution functions of phases are derived, and a number of theorems regarding general properties 
of phase distributions are proven. In Sec. IVII some of the analytic results are tested against numerical realizations of 
non-Gaussian fields. The complete form of the joint distribution function up to second order is given in Appendix lAl 



II. GAUSSIAN RANDOM FIELDS AND FOURIER MODES 

A. Fourier Transform and Power spectrum 

For a given random field f(x) in a <i-dimensional, Euclidean configuration space, x € M. d , the Fourier coefficient 
f(k) is defined by 

/>) = J d d xf{x)e- ik x . (1) 

Without loss of generality, we assume the field has a zero mean: (f(x)) — 0. The power spectrum, P(k), is defined 
by the relation 

(f (*)/(*')) = (2ir) d 6 d (k - k')P(k), (2) 

where S d is the c£-dimensional Dirac's delta function. The average (• • • ) represents an ensemble average. In an isotropic 
field, the power spectrum is an function of only an absolute length of the wavevector k. However, we do not assume 
the isotropy of space in principle, so that the configuration space can be anisotropic in general. The appearance of 
the Dirac's delta function in equation J2Jl is a consequence of statistical homogeneity of space. In this paper, we 
only assume the homogeneity of configuration space. Since the right hand side of equation (|2J) is non-zero only when 
k = fc', the power spectrum is always a real function. 

In this paper, we consider the field f(x) to be real: f(x) £ R. In this case, the Fourier coefficient of Eq. (JIJ satisfies 



f(fe) = /(-*)• 



(3) 
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Therefore, the power spectrum of Eq. J5J is equivalently given by 

(/»/>')) = (2n) d S d (k + k')P(k). (4) 

B. Probability Distribution Function of Fourier Modes in Random Gaussian Fields 

In random Gaussian fields, statistical properties of field distribution are completely specified by the two-point 
correlation function, 

£(*-*') = </(*)/(*')>. (5) 

In a homogeneous field, the two-point correlation function is a function of only a displacement of the two points. The 
power spectrum is a Fourier transform of the two-point correlation function (Wiener-Khinchin theorem): 

P(k) = J d d xe~ lkx £(x). (6) 

Therefore, statistical properties of random Gaussian fields are completely contained in the power spectrum, since 
the power spectrum and the two-point correlation function are mathematically equivalent. Higher-order moments in 
Gaussian fields are fully reduced to combinations of two-point correlation functions, and higher-order cumulants are 
zero: 



(f(x 1 )f(x 2 )---f(x N )) c = 0, (N>3). (7) 
The probability of taking particular values of a random Gaussian field f(x) is given by a functional, 



V[f(x)] = Aexp 



- 1 - J d&dtefwir^x-x'W) 



(8) 



where A is a normalization constant, and the function £ (x — x 1 ) is the "inverse" of the two-point correlation function, 
which is implicitly defined by 

/*•«.-. -x-v-o m 

Using the Fourier representation, this "inverse" correlation function is explicitly represented by 

r rfdh, ikx 

m-J&m (10 > 

The functional of Eq. (jHJ is a generalization of the multi-variate Gaussian distribution function to the continuum case. 
This functional contains all the statistical information of the field, and it is obvious from Eq. JSJ that the statistical 
properties of a random Gaussian field are fully described by the two-point correlation function, as claimed above. 

The distribution function of the Fourier coefficients is straightforwardly obtained from the Eq. (|SJ by changing the 
variables. Since the Fourier transform is a linear transform, the Jacobian of the transform is a constant. The integral 
in Eq. (jHJ) is given by 

A^/wr'c-^M^^M, (id 

where "uhs" indicates the integration over the independent modes taking the reality condition of Eq. © into account. 
Usually, one can take the "upper half sphere" , k z > of the fc-space. 

At this point, introducing a regularization box of volume V — L\L 2 ■ ■ ■ with periodic boundary condition is 
useful to discretize the Fourier space. The size of fundamental cells in Fourier space is given by Aki = 2-k/Li for 
the i-th direction (i = 1,2, ... ,d), and the right hand side of Eq. Ijlljl becomes a sum over discrete set of k in this 
representation. Defining the volume-normalized Fourier coefficient, 

h = W , (12) 
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the joint probability distribution function of the Fourier coefficients are given by 

uhs 

P[fh] =iIl eX P 
fe 

where A is a normalization constant, which differs from A by a Jacobian of the Fourier transform. The square bracket 
indicates the function of all Fourier modes in uhs, and not of only a particular mode. It is obvious in Eq. <|13|l that 
each Fourier mode in uhs is independent to each other in random Gaussian fields. Denoting /& = a& + ibk, where 
afc G K and bk € K are real and imaginary parts respectively, their distribution functions are given by 

P(a k ) = 1 



\.fk 



P(k) 



(13) 



a*; 



P(k) 



(14a) 



P(b k 



exp 



bk 

P(k)_ 



(14b) 



The real and imaginary parts of Fourier modes are independently distributed with rms P(k)/2 in Gaussian random 
fields. In terms of the polar representation, f k — |/fc|e i6,fc , the distribution function of each mode is given by 



P{\fk\,0 k )d\f k \d6 k = exp 



\h 



2 1 



P(k) 



2\f k \d\f k \ d6 k 
P(k) 2w 



(15) 



This distribution function does not depend on the Fourier phase 9 k , so that the Gaussian random fields always have 
random phases. The Fourier modulus \fk\ obeys the Rayleigh distributions in Gaussian random fields. 

Thus, the statistics of Fourier modes in Gaussian random fields are particularly simple: Distributions of every real 
and imaginary parts of Fourier modes are completely independent to each other, and their distributions are Gaussian. 
In other words, the distribution of Fourier modulus is given by a Rayleigh distribution and the Fourier phase is 
completely random. 

In what follows, we will investigate how this simplicity changes in non-Gaussian fields. We will see both independence 
of Fourier modes and Gaussianity of Fourier coefficients (and therefore randomness of phases) are broken in general. 



III. GENERAL DISTRIBUTION FUNCTION OF FOURIER MODES 

A. A formal expansion of non-Gaussian distribution functions 

One of the main purpose of this paper is to elucidate the general form of the joint probability distribution function 
of the Fourier coefficients P[fk\- For this purpose, we use the celebrated cumulant expansion theorem 0. In general, 
the cumulant expansion theorem states that the following identity holds for a random variable X: 

<e--)=ex P ff ^V> V (16) 

where (• • • ) indicates taking a usual average value, and (• • • ) c indicates taking cumulants. We put 

uhs 

X = Y,K k a k + L k b k , (17) 

k 

in the above equation, where Kk and Lk are just usual numbers, and variables a*, and bk are the real and imaginary 
parts of the Fourier coefficients, fk — ak + ibk («fc G M, 6fc e K). The summation is taken over any subset of 
independent Fourier modes. Since a mode with a wavevector k in a lower half sphere is always given by a mode with 
a wavevector —k in real fields, the wavevectors in the summation are taken only from uhs. 

With this substitution, the left hand side of Eq. IjlOII has the form of the Fourier transform of the probability 
distribution function P[a k ,b k ]. Therefore, performing inverse Fourier transforms, we obtain 

V[a kM = ff[ d ^^-e iX ^(±tf {xN) Y (18) 

J k \N=1 ' / 
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Since the original field fix) is assumed to have a zero mean, the first moments and cumulants vanishes: (a k ) = 
(ttfc)c = (bk) = (bk)c = 0- Thus the term of N = 1 in the last exponent of Eq. (|18|l vanishes. Substituting Eq. H17|) 
into Eq. 118J1 . using the binomial theorem, and expanding the sum over wavenumbers, we obtain 



V[a k ,bk] 



11 lis 

n 



dK k dL k 
2?r ~2tT 



x exp 



E 

m+n>2 



H) 



jhs uhs uhs 



In! 



'E^'EE'^'E^ ■■■a km bi 1 ■■■h n ) c K kl ■■■K km L h ■ ■ ■ U n 



k m l! 



In 



x exp 



uhs 



(19) 



where the sum over all non-negative integers m and n which satisfies m + n > 2 is taken. We single out the term of 
n + m = 2 in the first exponent, adopt substitutions K k — ► —id/da k , L k — * —id/db k in the rest terms, and perform 
the Fourier integrations. We obtain 



V[a k ,b k ] = exp 



E 

m+n>3 



(-1) 



m+n 



uhs uhs uhs 



uhs 



mm 



It,.! 



E " ■ E E " ■ E ^ afc i ■ ■ ■ ak m h ii ■ ■ ■ bi„ 



fei 



fem il 



d 



d d 
<9a fcl da km db kl 

where Pq is a Gaussian distribution function, given by products of Eq. I|14a|) and 114b|) : 



fe„ 



uhs _^ 



P(fc) 



V G [a k ,b h ], (2°) 



(21) 



The Eqs. (|20|l and (|21|l are the fundamental equations of expressing the non-Gaussian distribution function of Fourier 
coefficients in terms of higher-order cumulants. 

At this point, it is more convenient to use complex variables f k as an independent variables instead of a k and b k . 
One can consider simultaneous linear transformations of independent variables: 



h = a k +ib k (feeuhs) _ 
J-k = a k - ib k , 



(22) 



With this definitions, f k is defined in all k space, while a k and b k are defined in uhs space, so that the degrees of 
freedom is the same for both sets of variables. Carefully changing independent variables in the expression of Eq. I|2U|I . 
we obtain 



V\h\ = exp 



E 



both both 
,7V=3 " fci k N 



E---E •/<■• •••a- 



d 



d 



V G [fk], 



■V! 4- "" N ' c df kl df k , 

where the function P[f k ] and Pc[/fc] arc the formal probability distribution functions, defined by 



both 



uhs 



f[fh] \\ df k = V[a k , b k ] Y[ da k db kl 

k k 

both uhs 

Po[fk] \\ df k = Pc[ak, b k ] da k db k . 



(23) 



(24) 



(25) 



The summation in Eq. Il2.'{[| is taken over any subset of independent Fourier modes, provided that any two modes 
with ±fc are always included simultaneously ("both" stands for "uhs + lhs"). By calculating the Jacobian of the 
transformation of Eq. (|22|l . we obtain 



V[f k ] = 2- N ^V[a k ,b k ] 

uhs ^ 



f-kfk 



P{k) 



(26) 
(27) 
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The Eq. I|23|) with Eq. 1|27[1 can also be obtained by formally putting 



both 



X — ^2 Jkfk 



(28) 



in the Eq. (|16[) and formally following the similar calculations above as if fk's are real numbers. Therefore, the 
Eq. ill'MI) is considered as an analytic continuation of a general relation to the case of complex variables. In any case, 
the physical interpretation of the function P[fh] is given by Eq. (|24|l . the right hand side of which is a well-defined 
quantity. 

The higher-order cumulants in the exponent of Eq. 1|23|) are non-zero only when fci + • • ■ + k^ = because of the 
statistical homogeneity of the random field. In an infinitely large space, the polyspectra P^ N \ki, . . . , fejv-i), which 
are the higher-order counterparts of the power spectrum, are defined by cumulants of the Fourier coefficients as 



In finite- volume cases, 



f(ki)---f(k N ) 
= (2TT) d 6 d {k 1 + --- + k N )P ( - N ^k 1: . 



V l-N/ 2S K 



■+k, 



,kN-l) 



p( N \k u ...,k N ^), 



(29) 
(30) 



where 6 k K ^ is a Kronecker's delta defined by 



?(K) 



1, (fc = 0), 
0, (fe^O). 



(31) 



The polyspectra P( N '(ki, • • • , fcjv-i) are completely symmetric about permutations of arguments, and have symme- 
tries with respect to their arguments: 



pWt-fc!,... ,-fe w _ 1 ) = P w (fei,--- ,k N -i), 

pW(fei, fe 2 , • • • , k N ^) = P<< N \k 2 , • • • , fc w _i, -fci - k 2 fejv-i), 



(32) 
(33) 



and so on. 

Substituting Eq. J3UJ) into Eq. (23, we obtain 



V[f k ] = exp 



E 



/ i \N both both 



d 



N) V N / 2 



dfki df k „ 



Volfk] 



(34) 



This equation is a fundamental equation to derive the statistical distributions of Fourier coefficients in terms of 
higher-order polyspectra. One can formally take a limit of V — > oo in this expression, when all the Fourier modes are 
considered and included in summations. There are correspondences 



» {2n) d d d (k 1 
d d k 



(2tt)° 



dfk { ' 8f(k) 



■k N ), 



(35) 
(36) 

(37) 



where 6/8 f(k) is the functional derivative. In this limit, Eq. Ij34(l reduces to 



V[f(k)} =exp 



-J^T / d d k 1 ---d d k N (2n) d 8 d (k 1 + --- + k N )P^(k 1 , 

N=3 ' "* 



■ , few-i) 
5 

'6f(k 1 )'"~6f(k 



N 



r G [f(k)}, (38) 



7 



where P[f(k)] is a probability distribution functional of generally non-Gaussian random fields, and Pq [/(&)] is that 
of Gaussian random fields which share the same power spectrum with P[f(k)]. 
In the following, it will be convenient to introduce normalized variables, 



OLk 



fk 



which have a simple covariance matrix, 

{u k a k >) = <5fc+fc'- 

The probability distribution function of these normalized variables is given by 



V[a k ] = exp 



oo / 1 7y both both 

E^^E-E^---^) 



d 



d 



LN=3 



da kl da k „ 



V G [a k ] 



where 



P m (fci,..,M 



v /P(fc 1 )---P(fe A r_ 1 )P(fc JV ) 



are normalized polyspectra, and 



uhs ^ 

Pcbfe] =TT — cxp(-a_ fe a fc ) 



(39) 



(40) 



(41) 



(42) 



(43) 



The normalized polyspectra of Eq. I|42|l are non-zero only if k% -| fejv — 0, and satisfy the following relation 

= V 1 - N / 2 pW(k 1 ,...,k N ). (44) 



\a kl ■ ■ ■ a k „ 



The polyspectra do not depend explicitly on the volume V. This can be seen by the fact that the polyspectra are 
obtained by Fourier transforms of TV-point correlation functions: 



pW(k 1 ,...,k N -i) = Jd d Xl ---d d : 



XN-ie 



— i(fei -xiH hfcjv-i-aiiv-l 



where 



^ N) (x 1 - a; at, . . . , x N -i - x N ) = (f(xi) ■ ■ ■ f(x N -i)f(x N )) li 



(45) 



(46) 



is an iV-point correlation function. The case N = 2 of Eq. I|45|) is nothing but the old Wiener-Khinchin theorem 
of Eq. ©. Since JV-point correlation functions do not explicitly depend on the volume, the polyspectra also do not 
depend on the volume. 

Therefore, the expression of Eq. I|41|l can naturally be expanded by V^ 1 / 2 . After some algebra, we obtain 



V[a k ] 



oo oo 1 

l+ E y-«/ 2E _L E 



to! — (m+2)!---(n m + 2)! 

■■ ni,...,n m >l 



both 

E 



n±-\ \-rim=fi 

both 



^1 ni+2 1 n m +2 



where 



P, 



c9 



fci ■ ■ k n 



Po[ak}da kl da kl 



-V G [a k ] 



■P G [a k ], (47) 



(48) 
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H kl 


= a 


Hk t k 2 


= a 




= a 


Hkik^k^k^ 


= a 
+ 


kik^ksk&ks 


= a 
+ 


k2k 3 kiksks 


= a 



The last quantities are generalization of the Hermite polynomials 20, 21]. They are explicitly given by 

-fc 1 (49) 

-fextt-ha - <5fe 1+ fc 2 , (50) 

-fe x Q!_fe 2 a-fe 3 - [5f 1+k2 a- k3 + sym.(3)] , (51) 
-k x a-kvOi-k 3 a-k A - [^fc I +fc 2 a-fe 3 a-fe4 +sym.(6)] 

[^ 1+fe3 ^ 3+fe4 +sym.(3)], (52) 
_ fel a_ fe2 a_fe 3 a_fe 4 a_fe 5 - [5£ i+fe2 a_ fc3 a_ fc4 a_fc 5 +sym.(10)] 

K+* a ^+ fc4 tt-*.+sym-(15)], (53) 

Wfe... 

+ [Sk 1+ kjk 3 +k^-k^ ke + sym.(45)] - [<^ 1+fe2 <5£ 3+fe X 5 +fe 6 + sym.(15)] , (54) 

and so forth. In the above equations, "+ sym.(n)" indicates that n — 1 terms are added to symmetrize the preceding 
term with respect to fc's. For example, <5j^ i+fe2 a_fc 3 + sym.(3) = $% 1+k2 a-k 3 + S k i 2+k3 a^ kl + 8 k3+ki at-k 2 , etc - Tne 
Eq. (|47() is formally a series expansion by a dimensional quantity V~ x / 2 . The meaning of which is discussed in the 
following subsection. 

B. Explicit expansion up to second order 

Substituting the explicit representation of H klk2 ... into the expression of Eq. I|47|l . the general distribution function 
of the Fourier coefficients of a non-Gaussian field is obtained. In the following, we consider a Fourier amplitude A k 
and a Fourier phase 9k instead of a complex variable a k : 

a k = A *f\ ( fe e uhs ) . (55) 

where A k > and < 9k < 2tt. With these new variables, the distribution function of a Gaussian field is given by 
[cf. Eq. CS)] 

uhs uhs 

V G [A k ,9 k ]T[dA k d9 k =Y[2A k e- A <° 2 dAk— -. (56) 
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With these new variables, the Eq. (|47|l . up to second order in V 1 ^ 2 , reduces to 

_ uhs 

— = l+V- 1 ' 2 V A kl A k2 A k , cos (0 fcl + 9 k2 - 6 k3 ) (k u k 2 ,-k 3 ) 

{uhs _^ 
X] -Afe^^Afcj^cos^fe! +6» fc2 +6» fc3 - 6» fc4 )p (4) (fci,fc 2 ,fc3,-fe 4 ) 

uhs _^ 

+ ^2 -^ A k 1 A k2 A k ^A kl cofi(e kl +6 k . 2 -6 k3 -6 , fc 4 )P (4) [ki,k2,-k 3 ,-ki) 

k lt k 2 ,k 3 ,k4 
uhs 

- E 2 (^ 1 2 +A fc2 2 -l)p( 4 )(fc 1 ,fc 2 ,-fe 1 ,-fc 2 ) 

uhs ^ 

22 ^ A kt A k2 A k3 A kl A k5 A ke cos (9 kl +8 k2 ~ 8 k;i ) cos {9 kl + 8 k5 - 9 k& ) 



ki .k 2 



k 1 ,k 2 ,k 3 ,k4,k 5 ,k 6 



xj) (3) {k 1 ,k 2 ,~k 3 )p {3) {k 4 ,k 5 ,-k 6 ) 



ihs 



^ A fel A fe2 A fe3 A fe4 cos (0 fcl + 0fc 2 + 6> fc3 - fe J 

uhs 

X E^ 3) (^1, ^2, -fe 5 )p (3) (fe 3 , -^4, fc 5 ) 

uhs ^ 

-A kl A k2 A k3 A ki cos (0 fel + fc2 - fc; , - fe J 



k 1 ,k 2 ,k 3 ,k 4 



£ [p( 3 > (fci, fe 2 , -fe 5 )p (3) (fc 8) ki, -k 5 ) 

k 5 

+4p^ (fei, -fe 3l fc 5 )p (3) (fca, -*4, -fes) 



£ - (A kl 2 + A k2 2 + A k3 2 - 1) (fci.fca.-fca)] I + 0(y- 



3/2n 



(57) 



fel,fe 2 ,fe3 



where the wavevectors which appear as indices of and 0& are only summed over the uhs. 

The above expression, however, is still inconvenient for further investigations, because the modes in each summation 
with different fc-labels can be the same. We should expand each summation in Eq. I|57[) so that the different labels 
refer different modes. For example, in the second term the modes k\ and fc 2 can be the same and also can be different. 
In the same term, the modes k\ and k% should always different, because fe 2 = k% requires fc 2 = 0, which contradicts 
fc 2 G uhs. In the following we are going to integrate over some modes, in which case it is convenient when different 
labels refer to different modes. We carefully expand each summation in Eq. (|57() to summations in which the modes are 
mutually different. For example, the summation in the second term in Eq. (|57|1 should be separated into a summation 
with k\ = fc 2 and k\ ^ fc 2 , resulting in 

uhs 

^2 A kl A k2 A k3 cos (9 kl +8 k2 - 6» fc3 )p (3) (fci,fe 2 ,-fc 3 ) 

k%,k 2 ,k 3 

uhs 

= ^2 A k 2 A k2 cos{29 kl - 6» fc2 )p (3) (fci,fci,-fc 2 ) 

ki ,k 2 

uhs 

+ ^2' A kl A k2 A k3 cos (9 kl +9 k2 -9 k3 ) P (3 *> (k u k 2 ,-k 3 ), (58) 

ki,k 2 ,k 3 

where the summation Y^' k 8 ... indicates that all the modes are different to each other. On the left hand side, the 
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wavevector k% is identical to k\ + k 2 because of the Kronecker's delta in the normalized bispectrum. Therefore, the 
mode k 3 is automatically different from ki and k 2 because k%, k 2 ^ 0. The zero mode k = does not appear because 
it is a purely homogeneous mode which is excluded by imposing (f(x)) = 0. It is only when k\ — k 2 that two of three 
modes on the left hand side are identical. The case k\ = k 2 corresponds to the first term on the right hand side, and 
other case corresponds to the second term. 

Likewise, one can carefully expand other summations into summations of mutually different modes. After tedious, 
but straightforward manipulations, the result has a form, 

2 22 

W- = 1 + TV £ s ' (1) + 7 £ sf ' + ° {v "' % (59) 

v i— 1 i— 1 

where and are the summations given in Appendix lAl 



C. Meaning of the expansion by V ' 

We obtained the distribution function of the Fourier coefficients in a formal series by powers of V^ 1 ^ 2 . The volume 
V is a dimensional quantity, and the expansion seems meaningless if we naively choose the unit of length by V = 1 . 
However, powers of y -1 / 2 is always accompanied by the normalized polyspectra p^ n \ which are also dimensional 
quantities. The dimensionless quantities in our expansion is actually V 1_ ™/ 2 pW, which can be seen in Eq. I|41ll . 
Therefore, when these quantities are small for large n and higher-order terms are negligible, one can approximately 
obtain the distribution function by truncating the series, which is practically useful. 

In fact, the polyspectra are not dependent on the volume V, and so are the normalized polyspectra, when the 
wavevectors in the polyspectra are fixed. Thus the quantity V 1 ~ n ' 2 p^(ki, . . . , k n ) can be arbitrarily small by taking 
large V for a fixed set of wavevectors k%, . . . ,k n . Our expansion is actually an expansion by these quantities, assuming 
all the normalized polyspectra do not diverge. 

We should note that the number of possible wavevectors are large when the volume is large. This means that the 
number of terms which have a same order of V^ 1 ! 2 can be infinitely large if V — > oo. However, one does not have 
to consider all the possible modes in Eq. (|47|l . This equation is valid even in a case only a set of particular modes 
are considered, in which case the summation is taken over only wavevectors of that set of modes and the number of 
terms of a same order is finite. After the next section we consider distribution functions in which only particular set 
of modes are considered. 

Absolute values of the quantities V rl_ ™/ 2 p(™) are usually very small when the scales of wavevectors in the arguments 
of polyspectra are sufficiently smaller than the box size. We will see this property by some examples of non-Gaussian 
fields in section IVII below. Therefore, our expansion is actually efficient and turns out to be very useful for most of 
the cases. 



IV. TV-POINT DISTRIBUTION FUNCTIONS OF FOURIER MODES 



A. Zero-point Distributions 

As a consistency check, the probability distribution function P[Ak, 6k] in the form of Eq. I|59[) is shown to have a 
correct normalization, 

. uhs 



Jf[dA k d6 k -P = l, (60) 



even when the expansion is truncated. The integration on the right hand side can be explicitly performed using the 
expansion of V given in Appendix 17X1 with the Gaussian distribution Vq of Eq. I|56|) . The following integrals are useful: 



DC 



I(n) ee J A n ■ 2Ae- A dA = Y (l + - J , (61) 

\ — cos(n6 + a) = 0, (62) 
Jo 27r 

" g cosM + a) cosM + = { l mS(a ~ ^ [: J (63) 
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where n and m are non-negative integers. We use Eq. I|61|) with 1(0) = 1, 1(1) = 1(2) = 1, 1(3) = 3y/ir/4, 

1(A) = 2. Using these integrals, it is straightforward to show 



dA k d9 k Q { pV G = 0, 



for any mode k € uhs. The Eq. l|fiU)l is a direct consequence of this property. 



B. One-point Distributions 



(64) 



The distribution function of a particular Fourier mode f k can be obtained by integrating all modes but one in the 
general distribution function of Eq. (|59J) . The one-point distribution function Pi(A k ,9 k ) for a particular mode k is 
therefore given by 



r, uhs 

Pi(A k ,e k ) = / II dA k ,de k ,-v, 



(65) 



fc' 



where all the modes k' but a particular mode k are integrated over. 

Because of Eq. (I64|l . the Q± terms which have more than two different modes does not contribute to the above 
Eq. I|65|) . The only contribution comes from which has only one mode, and the Eq. (|65|l reduces to 



Pi(A k ,e k ) 



i 

V 



-A k 



Au 



?W (k,k,-k,-k) + 0(V- 2 ) 



P G (A k ,e k ) 



where 



P G (A k ,9 k )dA k d9 k = 2A k e- A * 2 dA k ^, 

Zn 



(66) 



(67) 



is a one-point distribution function of a Gaussian random field. 

There is a reason why only p^(k, k, — k, —k) contributes to the non-Gaussian correction and the terms with half- 
integer power of V do not contribute at all. The one-point distribution for a particular mode k is determined once 
the hierarchy of higher-order cumulants is provided, because of the cumulant expansion theorem. In the present case, 
we have two independent variables a k and ot- k for a single mode k £ uhs. Therefore all the cumulants have the form, 
(( a k) n («-fe) m ) c- Because of Eq. (|44ll . which is a manifestation of the translational invariance, such cumulants are 
non-zero only when n — m, and they are given by 



where 



((a k ) n (a_ k ) m ) c = 8 nm V x - n qV n \k), 
g( 2 ")(fc)=^)(fc,...,fc, -fc,...,-fc ), 

n elements n elements 



(68) 



(69) 



is a collapsed polyspectrum. Because the cumulants of Eq. (|68|l have integer powers of V the terms with half-integer 
power do not appear in the expansion of Eq. I|66() . 

In the case of one-point distribution function, obtaining higher-order terms in Ea. H66fl is not so difficult. To this 
end, we re-start the calculation from Eq. I|41fl . When we consider only a particular mode k 6 uhs, the sums over 
ki, . . . fejy in the exponent are only taken for k and — k. Because of the Kronecker's delta, the number of k and that 
of — k should be the same, so that N should be an even number, N — 2n. Taking proper combinatorial weights into 
account, we obtain a joint distribution function of a k and a^ k , 



P(a k ,a_ k ) 



1 

2^ 



exp 



E 



i 



,(2n) 



,n=2 v ' 



(k) 



d 2 



da k da- k 



exp(-a_ fc a fc ,) 



In terms of variables A k and 6 k defined in Eq. I|55|) , the differential operator in the exponent is given by 



d 2 



i / d 2 



1 d 



1 d 2 



da k da_ k 4 \dA 2 A k dA k A 2 d9 k < 



(70) 



(71) 
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The last factor in the new variables is given by exp(— A k 2 ) which does not depend on the phase 9 k . The Jacobian of 
the transformation is d(a k , / 'd(A k , 9 k ) = 2A k . Therefore, the general form of the one-point distribution function 
in terms of (A k , 9 k ) is given by 

Pi(A k ,9 k ) = — exp 

This equation provides a general expression of the one-point distribution function of a Fourier mode. There is an 
important theorem which is directly derived from Eq. I|72l) : 

Theorem 1 For a random field in a spatially homogeneous space, a one-point distribution function of Fourier phase 
P(9 k ) is always homogeneous: 

P(0k) = ^- (73) 



E 



7 (2") 



_^4( n !)2yn-l 



(k) 



( d 2 



\dA k 2 



_L d_ 

A k dA k 



exp(-A fc z 



(72) 



The proof of this theorem is trivial since the right hand side of Eq. i|72|) does not depend on the phase 9 k ■ The spatial 
homogeneity is crucial in this theorem. If the statistics of the random field is spatially inhomogeneous, the number of 
operator d/da k and that of d/dct-k does not necessarily agree because the lack of Kronecker's deltas in the exponent 
of Eq. (|41() . In which case, the exponent of Eq. I|72|l explicitly depends on the phase 9 k , and the theorem does not 
hold. 

The value of a phase 9 k determines positions of peaks and troughs of the Fourier mode fc, since the particular 
Fourier mode has the spatial dependence, f k e lkx — \ fk\e lk x+l6k , so that a spatial translation x — » Xq is equivalent 
to a phase shift 9 k — ► 9 k + k ■ Xq. Therefore, spatial homogeneity implies that there should not be preferred values in 
phases. This is an intuitive meaning of Theorem ^ However, when the other modes are simultaneously considered, 
the phases are not independently distributed and there are phase correlations among Fourier modes, as shown in the 
following sections. 

Another important consequence of Eq. JZIJ) is that the non-Gaussian corrections always vanish in the limit of V — > 
as long as all q<- 2n \k) are finite. Therefore, the following theorem holds: 

Theorem 2 For a random field in a spatially homogeneous space, a one-point distribution function of a Fourier mode 
approaches to be Gaussian when the spacial volume V is sufficiently large: 



P(\f k \,9 k )d\f k \d9 k ^exp 



provided that all the polyspectra of type P( 2n \k, 



l/fcl 



P(k) 



2\fk\d\f k \ d8 k 
P(k) 2ir ' 



when V — > oo, 



(74) 



-k) are finite for any positive integer n. 



This theorem is related to the central limit theorem. To illustrate the relation, we note the Fourier coefficients in the 
whole space are considered as superimposition of the Fourier coefficients of finite sub- volumes. A pair of sub- volumes 
which is sufficiently separated in space is almost independent to each other as long as the spatial correlations are not 
so strong on large scales. Expected separations of arbitrary pairs of such sub-volumes diverges in the limit V — > oo. 
The central limit theorem tells that superimposition of infinitely many independent random variables is normally 
distributed and the distribution function is Gaussian, irrespective to the distribution functions of the original random 
variables. The Fourier coefficients of the whole space is considered as superimposition of Fourier coefficients of many 
sub- volumes. Since these sub- volumes are almost independent, the distribution function of a Fourier coefficient of the 
whole space is expected to be Gaussian, even when the distribution in the sub-volumes are non-Gaussian. According 
to the Theorem |21 the one-point distribution function of a Fourier mode in a sufficiently large volume does not 
distinguish non-Gaussianity of the distribution. 

Expanding Eq. I|72() to arbitrary order in V~ x is straightforward. To second order, for example, we obtain 



Pi{A k ,9 k ) 
P G (A k ,9 k ) 



=1 



1 
1 

v 2 



1 



rA k 
1 

3G 
1 



A k 



A k e 



32 



A, 



-A k 



<z (4) (fc) 



W 

4 k 



(6) 



(k) 



3A k 2 



(k) \+0(V^). 



(75) 
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The first-order term agrees with Eq. I|66|) as it should be. It is straightforward to obtain higher-order terms by 
expanding Eq. (JJSJ. 

The distribution function of the original variables, fk = |/fc|exp(i#fc) is simply obtained from Eq. 1)75(1 . since the 
Jacobians of Pi and Pq are identical. In fact, Pi(|/fc|, Ok)/ PGi(\fk\, Ok) is just given by the right hand side of Eq. (|75|> 
with a substitution Ak = \fk\/ \/ P{k). 

C. Two-point Distributions 

The two-point distribution function for a couple of particular modes fci and k 2 is given by 

/uhs 
II dA k dO k ■ V, (76) 
k 

where all the modes fc but two particular modes fci and k 2 are integrated over. Instead of considering the full two-point 
distribution function of Eq. I|76[) . it is convenient to consider a reduced distribution function R2 defined by 

P2 (A kl ,9 kl ; A k2 , fe2 ) = Pi (A kl , 9 kl )Pi (A k2 , 9 k2 )[1 + R 2 {A kl , 6 kl ; A k2 ,6 k2 )]. (77) 

If there is no correlation between the two modes fci and &2, the reduced distribution function vanishes. The reduced 
distribution function represents the additional contribution that is not represented by one-point distribution functions. 
In the following, the reduced distribution function i?2(^4fci , kl \ Ak 2 , 0k 2 ) is abbreviated as R 2 (ki,k 2 ). The reduced 
distribution function does not change when the independent variables are changed since the Jacobian is the same for 
both sides of Eq. I|77|) . Therefore, the reduced distribution function for arbitrary set of independent variables, such 
as (Re/fc, Im/fc), (|/fe|, Ok), is simply obtained by represent the function in terms of the new variables. 

Most of the terms in Eq. I|59|l vanish in the Eq. (|76[1 . because of the integrations of Eqs. (|fci 1 1> — (|fci3|> . or Eq. I|64|) . The 
Qi terms with more than three different modes vanish in the integration of Eq. I|76[). The survived terms are only 
Q { i \ Q[ 2 \ Q { 2 2 \ Q { i\ and Qf\ The terms , Q {2) survive only when k 2 = 2k x or k x = 2k 2 , and the term Q {2) 

(2) (2) . (2) 

survives only when k 2 = 3fci or k\ = '3k 2 . The remaining terms Q\ , Q 2 always survive. However, the term Q\ 

does not contribute to the reduced function R 2 , because the contribution is absorbed in the one-point distribution 

(2) 

functions. In the term Q 2 , there is a symmetry k\ <-> k 2 in the summation, so that the symmetrization factor 
2 should be taken into account in the integration of Eq. I|76|) . As a result, the reduced function R 2 for two-point 
distributions is given by 

R 2 (k!,k 2 ) =-=A kl 2 A k2 cos(2<9 fcl - Ok.,)p (3) (fci, ki, -fe 2 ) + sym. (ki ^ k 2 ) 
W 



1 

2V 



A kl 4 A k2 2 cos 2 (20 kl -0 k2 )- 2A kl 2 A k2 2 - \a^ + 2A k 2 + A k 2 - 1 



(3) (fe 1 ,fe 1 ,-fe 2 ) +sym. (fei^fea) 
+ T^7^fci 3 ^fc 2 cos(36» fel - k2 )p {i) (fci, fci, fci, -k 2 ) + sym. (fci <-> fc 2 ) 
+ i {A k 2 Ak 2 2 - A k 2 - A k 2 + l)p (4) (fci.fe, -fci, -fc 2 ) 

+ 0{V- 3/2 ), (78) 

where ki,k 2 € uhs. The symbol sym. (fci <-> k 2 ) means an additional term which are needed to symmetrize each 
preceding term. The first and second terms of the lhs contribute only when k 2 — 2k\ or fci = 2fc2, and the third 
term contribute only when fci = 3fc2 or k 2 = 3fci. When the third term contributes, the first and second terms do 
not contribute, and vice versa. The forth term always contributes. 

The two-point distribution depends on phases only when the vector k 2 is proportional to fci. This property is 
understood by translational invariance of the statistical distribution. As considered in the previous subsection, phases 
appear in the distribution functions only in the form which is invariant under the phase shift Ok — » Ok + k ■ Xq. If k 2 
is not proportional to fci, there is not any way of making an invariant combination out of the two phases 0k ± and 0k 2 ■ 
If fc2 = cfci, a phase combination cOk x — 0k 2 is invariant. 
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D. Three-point distributions 

Next we consider the three-point distributions. The reduced three-point distribution function R 3 is defined by 

P 3 (fel,fe2,fe 3 ) = -Pl(fel) J Pl(fe2)i 3 l(fe3) 

x [1 + R 2 (k u k 2 ) + R 2 {k 2l k 3 ) + R 2 {k u k 3 ) + i? 3 (fc 1; fc 2 , fc 3 )] , (79) 

where we adopt notational abbreviations such as Pi (fci) = Pi(A] il , kl ), etc. Three modes, fci, k 2 , k 3 are all different 
from each other. The reduced distribution function is a component that is not represented by lower-point distribution 
functions. 

Contributions to the reduced three-point distribution function from Eq. l(59l come from , Q5 , Qg 2 \ 

( 2 \ ( 2 \ 

Q10 > Q\i • There are conditions among wavevectors under which these terms contribute to the distribution function. 

For example, the terms , Qg 2 ' survive only when there is a relation k 3 = k\ + k 2 or its permutation among the 
three modes, and so forth. The resulting reduced function R 3 is totally symmetric with respect to its arguments. It 

(a) 

is convenient to define an asymmetric function R 3 , with which the reduced function is obtained by 

i? 3 (fc l5 fc 2 , fc 3 ) = R<f\ki, k 2 , k 3 ) + R ( 3 \k 2 ,k 3 , fcx) + R { 3 \k 3 , fc 1? fc 2 ) 

+ Ri a) (k 2> fei, fc 3 ) + 4 a) ( fc i> fc 3, k 2 ) + R { 3 a) (fe 3 , k 2 , fci), (80) 

i.e., asymmetric functions with all the possible permutations are summed up to obtain the reduced function. The 
asymmetric functions can be extracted from expressions of Q\ given in the Appendix ^ resulting in 

4 a) (fc 1 ,fc 2 ,fc 3 ) = -^=A kl A k2 A k3 cos{9 kl +9 k2 -e k3 ) P {3) (fci.fe.-fcs) 
v V 

IpW (fci,fci,fc 2 ,-fe 3 ) 

,)p (4) (k u k u -k2,-k 3 ) 
1 



1 9 

—A k /A k2 A k3 cos (26 kl + k2 - 
1 „ 

^yA kl z A k2 A k3 cos (28 kl - 6 k2 



1 

V 



A fel 2 Afe 2 2 A fe , 2 cos 2 i 



- 9k 3 ) - \ {A kl 2 A k2 2 + A kl 2 A k3 2 + A k 2 A k 2 



~A kl — A k2 A k . 3 



^ (ki,k a ,-k 3 ) 



tj [A kl 2 A k2 3 A k3 cos (20 kl - 9 k2 )cos (26 k2 - 9 k3 ) - A k 2 A k2 A k . 3 cos (20 kl + 9 k2 - 8 k . 3 )] 



V 
1 

V 



x? (3) (fei,fei,-fe 2 )p (3) (k 2 ,k 2 ,-k 3 ) 
2A kl 3 A k2 2 A k3 cos (29 kl ~ 8 k2 )cos (0 kl + 8 k2 - 9 k:i ) - A kl 3 A k3 cos (30 fcl - 9 k . 3 ) 



2A kl A 2 k2 A k3 cos (9 kl - 29 k2 + 9 ka ) 



(k u k u ~k 2 ) P ^ (fcx.fea.-fca) 



+ (y- 3/2 ), 



(81) 



where ki,k 2 ,k 3 6 uhs. It is only when fci + k 2 = k 3 , or its permutations is satisfied that the first and fourth term 
on the rhs contribute. Similarly, 2ki + k 2 = k 3 for the second term to contribute, 2ki — k 2 + k 3 for the third 
term, 2ki — k 2 and 4fci = k 3 for the fifth term, 3fci = k 2 and 2ki = k 3 for the last term to contribute. The 
three wavevectors fci, k 2 , fc 3 should be in a 2-dimensional plane because there is at least one constraint among three 
wavevectors. In the last two terms, three wavevectors should be in a 1-dimensional line, i.e., three wavevectors should 
be parallel to each other. 

Specifically, if there is only a relation k\+k 2 = fc 3 , and if there is not any other relation among these wavenumbers, 
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the reduced three-point function is simply given by 

2 



R 3 (k 1: k 2 ,k 3 



A kl A k2 A k3 cos(9 kl +6 k2 -6»fe 3 )p (3) (k 1 ,k 2 ,-k. 



2 
V 



A kl z A k2 z A k / cos 2 {6 kl +6 k2 



0(V 3 ^ 2 ), (ki + k 2 = k 3 , no other independent relation). 



— - (A kl 2 A k2 2 + A kl 2 A k3 2 + A k2 2 A k3 ' 



-A k 2 - A k 2 - A k 2 + 1) 



p< 3 > (fci,fc 2 ,-fc 3 ) 



(82) 



The factor 2 relative to the asymmetric expression of Eq. (|S 1|) corresponds to a symmetry of fci k 2 . If there is 
another relation among the three wavevectors, the expression of Eq. I|82|l is insufficient and corresponding terms are 
added. 



E. Four-point distributions 

The four-point distributions can similarly be obtained. The reduced four-point distribution function R 4 is defined 

by 

P 4 {k l7 k 2l k 3 ,k 4 ) = P 1 {k 1 )P 1 {k 2 )P 1 {k 3 )P 1 {k 4 ) [1 + i? 2 (fc 1 ,fc 2 ) +sym.(6) 

+R 2 (k 1 ,k 2 )R 2 (k 3 ,k 4 )+sym.(3) + R 3 {k 1 ,k 2 ,k 3 )+sym.(A)+R 4 (k 1 ,k 2 ,k 3 ,k 4 )}, (83) 

where +sym.(n) indicates additional n — 1 terms which are necessary to symmetrize the preceding term with respect 
to ki, . . . ,ki. The same abbreviations for P 4 as in the case of P 3 , etc., are assumed. 

/ 2 \ ( 2 \ t 2 \ t 2 \ ( 2 \ 

Contributions to the reduced four-point distribution function from Eq. (|59l) come from Q 6 , Q7 > Q\ 3 \ Q\ 4 > 2i5 > 
Qi6 ) Qi7 ■ The term is absorbed in the non-reduced part of i? 2 • P-2 in Eq. Q83JI. As in the previous subsection, 
we define an asymmetric function , with which 

R 4 {ki,k 2 ,k 3 ,k 3 ) = 4 a) (fci,fc 2 ,fc 3 ,fe 4 )-l-sym.(24), (84) 
where the last term indicates additional 23 (= 4! — 1) terms to symmetrize the previous term with respect to fci, ... , k 4 . 
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The asymmetric function is given by 



R[ a) (k u k 2 ,k 3 ,k 4 ) = ^A kl Ak 2 A k3 Ak 4 cos(0 fcl +9 k2 +8 k3 -0 fc jp (4) (fci, fe 2 , fe 3 , -fe 4 ) 



1 

4V 
1 

V 



A kl A k2 Ak 3 A k4 cos(9 kl + 6 k2 - 9 k . 3 - 9 ki )p ( - 4) (fci,fc 2 , -fe3, -fe 4 ) 
2A kl 3 A k2 A k3 A ki cos (29 kl - hs ) cos (0 fcl - fc2 + fe J 



2A kl A k2 A k3 A ki cos (0 fel + 6» fe2 - fe;! - 9 ki 



,W (fci.fcx.-fcgJpW (fei,-fc 2 ,fc 4 ) 



1 

V 
1 

i 

r 



i 
r 



A kl 3 A k2 A k3 A k4 cos (20 fel - fe J cos (0 fel - fc2 - fe3 ) 

- A fel Afe 2 Afc3A fe4 cos(6» fel +<9 fc2 + fe3 -0 fe J p (3) (fei,fei,-fe 4 )p (3) (-fci,fc 2 ,fc 3 ) 
2A fcl 2 ^l fe2 A fe3 A fe4 2 cos (20 fel - 9 ki ) cos (0 fe2 - fe;! + fc J 

- A fcl 2 v4 fc2 Afc 3 cos (20 fel + fc2 - fc3 )lp (3) (fci, fci, -fc 4 ) p (3) (fc 2 , -fc 3 , fe 4 ) 



A fcl 2 A fc2 A fe3 Afc 4 2 cos (20 fel - 0fc 4 ) cos (0fc 2 + fc 
1 



- -A fel 2 A fc2 A fc3 cos (20 fel - fc2 - fe J 



(3) 



(fci,fci,-fc 4 )p (3) (k 2 ,k 3 , -fe 4 ) 



4A fcl 2 A fc2 A fc3 A fc4 2 cos (0 fel + fc2 - fc4 ) cos (9 kl - 9 k3 + 9 ki ) 



— 2A kl 2 A k2 A k3 cos (20fej + fc2 — 9 k . 3 ) — 2A k2 A k3 A ki 2 cos (0 fe2 + 9 k3 — 29 ki ) 



x p 



(3) 



(fci,fc 2 ,-fc 4 )p (3) (fei,-fe 3 ,fc 4 ) 



(85) 



where fci, fe 2 , fe 3 , fe 4 € uhs. 

Specifically, if there is only a relation fci + fc 2 + k 3 = fe 4 , and if there is not any other relation among these 
wavevectors, the reduced four-point function is simply given by 



J? 4 (fci,fe 2 ,fc3,fe4) = —A kl A k2 A k3 A ki cos(0 fel + 0fc 2 + 0fe 3 - 0fc 4 )p (4) (fci,fc 2 ,fc3, -fc 4 ) 
+ 0(y~ 3 ^ 2 ), (fci + fc 2 + k 3 = fc 4 , no other independent relation). 



(86) 



The symmetry factor 6 of permutating (fci,fc 2 ,fc3) is taken into account. When there is an additional relation, such 
as 2fci = fc 4 etc., a corresponding term, such as the fourth term in Eq. (f%5)) . should be added to the above equation. 
Similarly, if there is only a relation fci +fc 2 = fc3 + fc 4 , and if there is not any other relation among these wavenumbers, 
the reduced four-point function is given by 



i? 4 (fci,fc 2 ,fc 3 ,fc 4 ) = ^-A kl A k2 A k3 A kl cos(0 fcl + fc2 - fc3 - 0fc 4 )p (4) (fci,fc 2 ,-fc 3 , -fc 



V 



0(V 3 / 2 ), (fci + fc 2 = fc3 + fc 4 , no other independent relation). 



(87) 



F. Higher-point distributions 

It is also straightforward to obtain expressions of five- and six-point distribution functions from Eqs. (| A2 If) — A25p . 
Up to the order V -1 , five- and six-point reduced distributions are present only when the wavevectors satisfy at least 
two conditions. For example, the five modes fci, . . . , fcs with a condition fci + fc 2 + fc3 + fc 4 = fcs do not have a reduced 
five-point distribution function up to this order. This is consistent with the fact that iV-point correlations of Fourier 
modes are of order I/ 1-JV / 2 as given by Eq. (|30|l . 
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V. PHASE CORRELATIONS 



One of the most prominent feature in non-Gaussian fields is that the Fourier phases are not completely random. 
Therefore, the non-Gaussianity is sometimes called as "phase correlations". However, explicit forms of the phase 
correlations in non-Gaussian fields was unknown. Unlike in the random Gaussian fields, Fourier moduli and Fourier 
phases are not independently distributed, as seen from the distribution function we derived above. Thus, correlations 
only among phases does not have enough information on non-Gaussianity. Moduli and phases of different modes are 
mutually correlated. 

When one is interested only in phase distributions, the distribution function of phases is obtained by integrating over 
Fourier moduli. Such integrations are simply performed by using Eq. 1(61(1 . In the expressions for Qj in Appendix lAl 
the integrations result in replacements 

V-l, V-^, V-2, (88) 

since all the different labels of wavevectors corresponds to actually different modes. The same replacements are applied 
in the expressions of the ./V-point distribution functions. The one-point phase distribution function is constant because 
of Theorem ^ 

2nP 1 (9 k ) = 1. (89) 

This can be explicitly confirmed by Eq. ((66(1 and also by Eq. (|75|l . 

For higher-point distribution functions of phases, reduced functions are defined along the lines of defining Rn- From 
the iV-point distribution function of phases Pjv(^fen ■ • ■ the reduced distribution functions Cn are iteratively 

defined by 

{2n) 2 p 2 {e kl ,e k2 ) =i + c 2 (e kl ,e k2 ) (90) 

(2n) 3 P 3 (9 kl , 9 k2 , 9 k3 ) =1 + C 2 {9 kl ,6 k2 ) + C 2 (0 k2 ,9 ks ) + C 2 (6 k3 ,0 kl ) + C 3 (8 kl , 8 k2 , 8 k3 ) (91) 
(2tt) 4 P 4 (0 fcl , • ■ • , ki ) =1 + C 2 {8 kl ,9 k2 ) + sym. (6) + C 2 (9 kl , 9 k2 )C 2 {9 ks ,9 ki ) + sym. (3) 

+ C 3 {9 kl , 6 k2 ,9 k3 )+ sym. (4) + C 4 (6 hl , 6 k2 , 9 ka , fe J (92) 

Such reduced functions are just given by integrals of the function R^: 

C N (6 kl ,...,9 kN ) = (2n) N J dA kl ---dA kN P 1 (ki)---Pi(k N )R N (k 1 ,...,k N ). (93) 



A. Two-point distributions 



The two-point distribution functions of phases are obtained from Eqs. 1)66(1 . 1(78(1 and ((93(1 . Up to the order U" 1 , 
the one-point function Pi(k) in Eq. 1(93(1 can be replaced by a Gaussian function Ps(k). Thus the reduced function 
C 2 is just given by replacements of Eq. (|8*5j) in Eq. lf78)l. As a result, the fourth term of Eq. l(T5l) vanishes and the 
phase correlations between two modes appear only when one of the modes is parallel to the other and the proportional 
factor is either 2 or 3. Specifically, when 2ki = k 2 = k, 



C 2 



Ik, 2k 



cos(2fl fc - 9 2k )p^ (fe, k, -2k) + -- cos [2(26 k - 9 2k )} p (A > (fc, k, -2k 



2y/V 

and when 3fci = k 2 = k, 



(3) 



2U 



C 2 (0 k ,03k) = cos(30 fc - 9 3k ) P ^ (fc, k, k, -3k) + 0{V-' 6/ ' z ). 
o V 



-3/2-1 



0(V~ 3/2 ), (94) 



(95) 



Otherwise, there is not any other two-point correlation between phases up to this order. In the above two cases, 
the phase correlations are present only between two modes in which wavevectors are mutually parallel and their 
proportional factor is a simple integer. This property is not specific for the lower-order expansion. In fact, there is a 
following theorem: 
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Theorem 3 Phase correlations between two different modes with wavevectors k± and h,2, where k\,k,2 € uhs, are 
present only when the two wavevectors are parallel to each other, and the proportional factor c is a rational number. 

Proof. The two-point distribution function is obtained from Eq. I|47|l where the wavevectors in the summation take 
their values only from four vectors, fei, — k%, k2, and — k2. Let the number of these wavevectors in a term be m, m\, 
n 2 , m2, respectively. Because of Kronecker's delta's in the polyspectra, the equation n\k\ — m\k\ + ri2fc2 — m^k2 = 
should be satisfied in the term. Therefore, we obtain 

m - mi , . 

k 2 = ki (96) 

m 2 — n 2 

unless n\ = mi and n 2 = ni2 are simultaneously satisfied. The theorem is proven if a term in which the conditions 
n\ = mi and n2 — m2 are satisfied does not contain phase distributions. In the term with n\ — m\ and n 2 = 1112, 

the generalized Hermite polynomial has a form of H kl kl ... k2 k 2 , where the numbers of ki and —ki are the same, 

and that of k2 and — fe 2 are the same. In explicit forms of the polynomials, as illustrated in Eqs. (|49|) - (|54|l . the 
terms with Kronecker's delta's, S£ +fel , S K ki _ ki , <5j? 2+fc2 , $ K k2 - k2 , <5]? 1+fc2 , £* fel _ fe2 , are zero since fei,fe 2 G uhs. Only 
terms with _ k and S_ k , fea survive. Therefore, the numbers of fci and — k\ are the same and that of fc 2 and — &2 
are the same even in each term of the generalized Hermite polynomials. That means each term is a function of only 
combinations of a kl a>- kl = \ct kl \ 2 and a k2 a>- k2 — \ak 2 \ 2 , which do not depend on phases. ■ 
According to this theorem, there are phase correlations between modes fci and fc2 only when a relation nk± = mk2 
is satisfied for some different integers n, m. In terms of music, this is a relation of "overtones" . Considering a 
fundamental tone k, the two modes are given by m-th and n-th overtones, k\ = mk, k2 = nk. In the Eq. Ij94|l . 
for example, these numbers are m = 1 and n = 2, and they correspond to a fundamental tone and a first overtone 
(an octave), respectively. In the Eq. I|95|l . the relation is a perfect twelfth (an octave and fifth). A phase correlation 
between 2k and 3fc corresponds to "a perfect fifth" , and that between 3fc and 4fc corresponds to "a perfect fourth" , 
and so on. 

In Eqs. (|94(l and Ij95(l . lowest-order deviations from random distribution of phases are of order C(y -1 / 2 ) for a 
two-point phase distribution of modes k and 2k, and ©(U -1 ) for modes k and 3fc. The deviations are of higher order 
when a proportional factor of the two wavevectors is not simple. More precisely, we have a following theorem: 

Theorem 4 Deviations from a random distribution of phases for a two-point phase distribution of modes k\ and k2 
with a relation nki = m,k 2 , where n and m are irreducible positive integers, are of order 0(V~( n+rn ~ 2 ^ 2 ), or higher. 



Proof. In the proof of Theorem |31 the term specified by positive numbers m, mi, n2, 7712 has a order of 
Q(y-{n 1 +m 1 +n 2 +Tn 2 -2rn')/2^ w ^ eIe m ' [ s ^ numD er of polyspectra in the term in Eq. g7|). Firstly, the lowest- 
order contribution comes from the terms which has only one polyspectrum, m' = 1. Secondly, from Eq. (|96(l in 
the context of the proof of Theorem |3 n/m — (ni — m\)j(m2 — ^2) should be satisfied to have non-trivial phase 
distributions. Under this constraint, the lowest-order term is given by (ni, mi, n2, TO2) = {n, 0, 0, m), (0, n, m, 0), and 
the order is 0(V( 2 ~ n ~ m ^ 2 ). This term is accompanied by a specific polyspectrum. If the corresponding polyspectrum 
vanishes, the lowest order is higher than that. ■ 
In fact, the lowest-order term in this theorem can be explicitly obtained by inspecting Eq. H47|) . From the proof 
above, and counting overlapping factors in Eq. (|47l) . we obtain 



y(2-n-m)/2 
7 1 m I 



C 2 (0 kl ,0k 2 ) = - , r- P ( - n+m) (k 1 ,...,k 1 ,~k2,...,-k 2 ) 



x J 2A kl dA kl 2Ak 2 dAk 2 (H kl ,. . .,k u -k 2 , ...,-h 3 + H- kl , . . . . -k t ,k 2 , ■ ■ -,k 2 ) CX P {-Ak ± 2 - A k 

Since k\ and &2 are in uhs and are mutually different, all the Kronecker's delta's are zero in expansions of generalized 
Hermite polynomials appearing in the above equation. Therefore, we have 

Hk lt ... , fei,-fc 2 , ■■ -k 2 + -ff-fci , -k x ,k 2 , ...,k 2 — (a-fci)"(afe 2 ) m + (a.k 1 ) n ( y ot-k 2 ) m (98) 

= 2 A kl n A k2 m cos (n6 kl - m6 k2 ), (99) 
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and Eq. (|9"T|) reduces to 

OT/(2-n-m)/2 / n\ / m\ 

c 2 (e kl ,e k2 ) = — — — r (1 + -J r (1 + - J cos ( n e kl - m e k2 ) 

x p< n+m >(fei, . . . , fei, -fe, . . . , -fe) + 0(y( 1 -"-™)/ 2 ), 
> „ * „ ' 

n m 

(nki = mk 2 for irreducible integers m,n). (100) 

The lowest-order correction terms in Eqs. (|94|l and (|95|l are reproduced from Eq. (|10L)f> by putting m = l,n = 2 
and to = l,n — 3 respectively. Lowest-order corrections of all the other two-point phase correlations are generally 
obtained by Eq. (jTUU|) . 

In terms of music again, the phase correlations are stronger for modes with a simpler overtone relation. For 
example, a pair of modes with an octave, (m,n) = (1,2), results in a phase correlation of order V~ x ^ 2 . That 
with a perfect twelfth (an octave and fifth), (m,n) = (1,3) results in a correlation of order V -1 , Similarly, double 
octaves, (m,n) — (1,4), and a perfect fifth, (m,n) = (2,3) result in correlations of order V~ 3 ^ 2 . A perfect fourth, 
(m,n) = (3,4), corresponds to the order V~ 5 ^ 2 , and so forth. In other words, Fourier phases are correlated between 
modes when the modes have an harmony. The simpler the harmony sounds, the stronger the phase correlation is. 



B. Three-point distributions 

The three-point distribution functions of phases are obtained from Eqs. I |66| l. Q8U[). I|81|l and l|93|) . Up to the order 
V -1 , the reduced function C 3 is just obtained by replacements of Eq. I|88() in Eq. i|81|l and symmetrization. The 
phase correlations among three modes appear only when the three wavevectors are linearly related. Specifically, when 
k\ + &2 = &3 and there is not any other relation among these wavevectors, integrations of Eq. (|82|) result in 



C 3 (fci, k 2 , k 3 



7T 3 / 2 



cos((9 fcl +9 k2 -6 k3 )p {3) (k u k 2 ,-k a 



+ -cos[2(9 kl +0 k2 -6 k3 )} |pW (ki,k 2 ,-k 3 ) 
( ki + k 2 = &3, no other independent relation). 



<D(V 



-3/2n 



Similarly, it is easy to obtain phase correlations when there is only one linear relation of either 2fei + k 2 
2ki = k 2 + k 3 . For the first case, 



(101) 

k 3 or 



C 3 (fei,fe 2 ,fe 3 ) = — cos(20 fcl +6 k2 -9 k3 )p^ (k 1 ,k 1 ,k 2 ,-k 3 ) + 0{V 
(2ki + k 2 = k 3 , no other independent relation), 



-3/2N 



(102) 



and for the second case, 



C 3 (kx,k 2 ,k 3 ) = ^cos(28 kl -6 k2 -0 fc3 )P (4) (fei.fei, ~k 2 -k 3 )+0{V~ 3 > 2 ) 
(2fci = k 2 + k 3 , no other independent relation). 



(103) 



Another specific example is that the three modes are given by k, 2k, 3k. In this case, the first, third, fourth, and 
fifth terms in symmetrized Eq. (|81|l are relevant. The phase correlation in this case is calculated to be 



C 3 (fe,2fc,3fe) 



7T 3 / 2 

1 



cos (0 fc + e 2k - e 3k ) P ^ (fe, 2k, -3k) 



cos [2 (9 k + 6 2k - 9 3k )} (k, 2k, -3k) 
cos (8 k - 29 2k + 9 3k ) fopW (2k, 2k, -k, -3k) 

8 V 



-p (3) (k, k, -2k)p (3) (k, 2k, -3k) 



0(V 



-3/2n 



(104) 
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To obtain the three-point distribution P 3 (0k, Oik-, 3 k) from the above expression, it should be noted that the two-point 
contributions given by Eq. I|94|) and (|95|l have to be taken into account. 

Yet another example is the case that the three modes are given by k, 2k, 4k. In this case, 

C 3 (fc, 2k, 4k) = ^- cos (28 k + 2k - 4 fe)i> (4) (k, k, 2k, -4k) 

+ ^- [3 cos (20 k - 2k ) cos (202k - 04k) - 2 cos (20 k + 2k - 4k )} 
o V 

x p ( V (fc, k, -2k)p (3) (2k, 2k, -4k) 
+ 0{V~ 3/2 ). (105) 

The Theorems [3]and0]can be generalized to the case of three-point distributions. Corresponding to the Theorem^ 
we have a following theorem: 

Theorem 5 Phase correlations among three different modes with wavevectors k±, k2 and k 3 , where k\,k2,k 3 € uhs ; 
can be present only when there is a linear relation with integral coefficients among the three wavevectors: 

jik! + j 2 k 2 + j 3 k 3 = 0, (106) 

where ji (i = 1,2,3) are integers and at least two of ' j±, j 2 , js are not zero. 

Proof. The proof of this theorem is quite similar to that of Theorem [3J We consider the Eq. I|47|l where there are 
six vectors, ±fei, ±fe2, ±£3. These wavevectors should satisfy n\ki — m\k\ + n 2 k 2 — m 2 k 2 + ^3^3 — m 3 k 3 = in 
each term, where n\ is a number of k\, ui\ is a number of —k\, etc., in the corresponding term. Therefore, Eq. I|10tj|) 
should be satisfied with ji = rii — m,-. If only one of jVs is not zero and the other two are zero, the corresponding 
wavevector is identically zero, which is excluded by the assumption. Therefore, the theorem is proven if the terms 
with ji = j 2 = js = do not have phase dependence. The rest of the proof is almost the same as that of Theorem [3] 
■ 

It is immediately follows from this theorem that the three wavevectors are laid on a two-dimensional plane even 
when the spatial dimension is more than two. In other words, if the three wavevectors are not laid on a common 
plane, the phases are not correlated at all. 

When there is only one linear relation among the three wavevectors, as in the cases of Eq. <|101[) . I|102fl . (|103fl . 
the lowest-order term of C 3 has a simple form. In fact, this simplicity is a general property, which can be shown 
by generalizing the derivation of Eq. (|100fl to the case of three-point distributions. We consider the case that the 
only linear relation is given by Eq. (|106f) where j\, j 2 , j 3 are non-zero, irreducible integers. We define the signs of ji 
(i = 1,2, 3) by 

Although there is an ambiguity of choosing an overall sign in Eq. (|106|) , the following arguments are not affected by 
the choice. Counting overlapping factors in Eq. 147(1 . the lowest-order contributions are given by similar calculation 
of Eqs. (|97 |) ~ (|TUUI) : 

^,e kM - Mmm r ^1 + — j 1^1 + - j r + -j 

x cos(ji0 kl + ]20k 2 + h0k 3 )p (bl ^ +b ^ +bM {siki, . . .,s 1 k 1 ,s 2 k 2 , . . .,s 2 k 2 ,s 3 k 3 , . . .,s 3 k 3 ) 

" v ' " v ' v ' 

bil U2I Us I 

+ on/(Mjii-U2i-U3i)/2\ 

O'i^i + jiki +i3^3 = 0, jij2j3 7^ 0, no other independent relation). (108) 

The "unconnected part" of Eq. (|91fl . i.e., contributions of C2 on the left hand side, are not present in this case, because 
it is assumed that there is not any linear relation between two of the three wavevectors. It is crucial that any of j\, 
j 2 j 3 is assumed to be non-zero in the above equation. 

If the spatial dimension is one, the number of independent linear relation among three wavevectors can not be only 
one, and the Eq. 1)1081) does not apply at all. For example, Eq. (|105H is applicable even in one dimensional space, in 
which case the lowest-order term is not given by just one polyspectrum. 

If the spatial dimension is two, any three wavevectors can not be linearly independent. There is always a linear 
relation among three vectors. When any two of the three vectors are not parallel to each other, there is only one 
linear relation. The simpler the coefficients of the linear relation is, the larger the three-point phase correlation is. 
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The lowest-order terms in Eqs. ()101|l - (|103(l are reproduced by putting (ji,j2,js) = (1, 1, — 1), (2,1,-1) and 
(2,-1,-1), respectively. There are more than one linear relation in Eqs. (|1L)4|) and l|l(J5|) . and the Eq. (|99|) is 
not applied to them. 



C. Four-point distributions 



The four-point distribution functions of phases are similarly analyzed as the three-point functions. The phase 
correlations among four modes appear only when the four wavevectors are linearly related. Specifically, when fci + 
k 2 + k% = fc 4 and there is not any other relation among these wavevectors, integrations of Eq. (|86() result in 

_2 

C4{e kl ,e k2 ,e k3 ,6 k4 ) = —cos(8 kl +e k2 +e k3 -e k4 )p (4) (fci,fc 2 ,fc 3 ,-fc 4 ) + o(v~ 3 / 2 ), 

8 V 

(fci + k 2 + ^3 = fe 4l no other independent relation). (109) 

When ki + k2 — + fc 4 and there is not any other relation among these wavevectors, integrations of Eq. 1(87(1 result 
in 

c A (e kl ,e k2 ,e k3 ,e k j = ^cos(e kl +e k2 -e k3 - 9 ki ) P ^ (k^k^-k^-k^ + o(v- 3 / 2 ), 

(fci + &2 = fc 3 + fc 4 , no other independent relation). (110) 

When there are more than one linear relation, the phase correlations have different expressions. The explicit 
expressions are derived from integrations of the Eq. (|85fl and symmetrization. The seventh term in Eq. I|85|) identically 
vanishes by integrating over the moduli A kl , . . . ,A k4 , and does not contribute to the phase correlations. 

For example, we consider a case that 2fc 2 — k 3 — and fci — k 2 + fc 4 = are simultaneously satisfied. These 
equations are equivalent to the expressions fc 3 = 2fci and fc 4 = fc 2 — fci. The third term in Eq. 1)85(1 contributes to the 
phase correlations in this case. Additionally, an equation fci + fc 2 — fc 3 — fc 4 = is also satisfied, and the second term 
contributes as well. Other terms and their symmetrization are incompatible with the conditions. Therefore, we have 

C4(9 kl ,9 k2 ,9 k3 ,e k j = — cos(e kl +e k2 -e k3 -e ki )p {4) (fc^fo.-fca,-**) 

ov 

+ ^7 I 3 cos ( 26l fei - 0fc 3 ) cos (9 kl 9 k2 + 8 ki ) - 2 cos (6 kl + 9 k . 2 - 8 ks - 9 ki )] 

x (fci, k lt -fe 3 )p (3) (fci, -fe 2 , k 4 ) + 0(V~ 3 / 2 ), 
(2fei = fe 3 and ki + fc 4 = k 2 , no other independent relation). (HI) 

Similarly, considering other cases when either fourth, fifth, or sixth term in Eq. (|85() is effective, we obtain 

7T 2 

Ci{9 kl ,9 k2 ,9 k3 ,8 ki ) = —cos(8 kl + 9 k2 - 9 ks -8 ki )p {A) (k 1 , k 2 , -fe 3 , -fe 4 ) 
8 V 



n 2 



[3 cos (29 kl - 9 ki ) cos (6 kl - 9 k2 - 9 ks ) - 2 cos {6 kl + 9 k2 + 9 ks - 9 ki )] 



16V 

x p& (k u kx, -fc 4 ) P {3) (-fci, fca, fcs) + 0(V- 3 / 2 ), 
(2fei = fc 4 and fci = fc 2 + fc 3 , no other independent relation), (112) 

C A (9 kl ,9 k2 ,9 k3 ,e ki ) = ^cos(29 kl -9 k2 +9 k3 -29 k Jp^ (fc 1; fc 1; -fc 4 )p (3) (fc 2 ,-fc 3 ,fc 4 ) 
+ C(^- 3/2 ), 

(2fci = fc 4 and fc 2 + fc 4 = fc 3 , no other independent relation), (113) 
C 4 (^ fel , 9 k2 ,9 k3 ,6 ki ) = ^ cos (20 fcl + 8 k2 + 9 ka - 20 fc Jp (3) (fc x , fc 1} -fc 4 )p^ (fc 2 , fc 3 , -fc 4 ) 

+ 0(y- 3 / 2 ), 

(2fci = fc 4 and fc 2 + fc 3 = fc 4 , no other independent relation). (114) 
There is a theorem which is similar to Theorem If there is only one linear relation among the four wavevectors, 
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the lowest-order term of C4 can also be obtained as a trivial generalization of Eq. i|108|l : 

2y( 2 -|j'il-b'2|-|33|-|j4l)/2 



|ji|!|j 2 |!|j 3 |!|j4| 



X r(i + M)r(i + M)r(i + M)r(i + M 

x cos(ji6» fcl +j 2 k2 +j3&k 3 + jiO ki ) 

x p(b'il+b'=l+b3|+b4|)( Slfelj ^ sifei, s 2 fe2, ■ ■ ■ , s 2 fei, s 3 fe3, ■ ■ ■ , s 3 fci, S4J24, ■ ■ ■ , 

" v ' " * ' " v ' " v ' 

b'i| lis I lis I b'4| 

_l_ 0(y(i-|id-li2l-li3|-li4|)/2\ 
(jifei + j 2 fc 2 +^3^3 +J4fe4 = 0, jijzjsjA, 7^ 0, no other independent relation). (115) 

Eqs. H109fl and (|110|l are reproduced from the above formula. 

If the spatial dimension is two or smaller, the number of independent linear relations among four vectors can not 
be just one, and the Eq. I|115l) does not apply. For example, Eqs. Hlll|) ~ (|114l) is applicable even in two dimensions, 
and the lowest-order term is not given by just one polyspectrum. 

If the spatial dimension is three, any four vectors cannot be linearly independent. There is always a linear relation 
among four vectors. When any three of four vectors do not have a common plane, there is only one linear relation. 
The simpler the coefficients of the linear relation is, the larger the four-point phase correlation is. 



D. TV-point distributions 

Now it is straightforward to generalize the Theorem [3] and Eqs. I|99|l and (|115|l to cases of TV-point distributions of 
phases, where TV > 5. As a counterpart of Theorem^ we have the following theorem: 

Theorem 6 Phase correlations among TV different modes with wavevectors ki, fc 2 , • ■ • , fc/Vj where TV > 3 and 
fei, fe 2 , . . . , fejy € uhs, can be present only when there is a linear relation with integral coefficients among the TV 
wavevectors: 

jifel + ]2k 2 H h ]Nk N = 0, (116) 

where ji (i = 1, 2, . . . ,N) are integers and at least two of j\, j%, . . . , Jn are not zero. 
The generalization of Eqs. I|108(l and l|115|) to TV-point distributions is given by 

2y(2-lii|-li2| li«|)/2 



Cn (#fci , 0k 2 > ■ ■ ■ j ®k N ) — 



lii|!|i 2 |!---|iiv|! 



x r (i + r (i + • • • r (1 + cos { n e kl + j 2 e k2 + ■■■+ j N e kN ) 

x p(bd+b'2l+-+b>l)( Sl fe l! . 5 Sl fe lj S2 k 2 , . . . , S2 fei, . . . , s N k N , s N k N ) 

S v ' S v ' S * ' 

lii I \h I Ij'n I 

+ ^(^^"bd-lial \3n\)/2\ 

(jifci +J 2 fe 2 + • • • + ijvfejv = Oj Jij2 • • • Jjv 7^ 0, no other independent relation). (H?) 

However, if the spatial dimension d is TV — 2 or smaller, the number of independent linear relations cannot be just 
one, and the above equation does not apply for TV > d + 2. In three dimensional space, TV > 5 of Eq. I|117[l does not 
apply. 



VI. NUMERICAL EXAMPLES 



In this section, analytic results derived in previous sections are compared with numerical realizations of simple 
non-Gaussian fields. One of the purposes of this section is to check the analytic results, which is derived above by 
lengthy calculations. Another purpose is to show how accurately truncated expressions in series of V^ 1 !" 1 reproduce 
the numerical distributions of Fourier modes. 
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FIG. 1: An example of a Voronoi tessellation in 2-D space. 



We artificially generate realizations of non-Gaussian fields on regular sites in a three-dimensional cubic box. The 
number of sites, which corresponds to the resolution of the fields, are 256 3 . We generate three types of non-Gaussian 
fields, i.e., Voronoi tessellation fields, lognormal fields, and quadratic Gaussian fields. The first one is purely non- 
Gaussian, and the last two fields are generated by simple non-linear transformations of Gaussian fields. 



Our first examples of non-Gaussian fields are Voronoi tessellation fields. To define a concept of Voronoi tessellation 
fields, first we consider a set of points S which are randomly distributed in space. A Voronoi tessellation [2^ is a 
set of cells, each of which delimits a part of space that is closer to a point of S than any other point of S. An 
example of a Voronoi tessellation in 2-dimensional space is illustrated in Fig. ^ In 2-dimensional space, a Voronoi 
tessellation consists of 1-dimensional lines on which the number of nearest points of iS is more than one. Similarly, in 
d-dimensional space, it consists of (d — l)-dimensional hypersurfaces on which the number of nearest points of S is 
more than one. In 3-dimensional space, it consists of 2-dimensional surfaces, which we call Voronoi surfaces below. 

A Voronoi tessellation itself is not a random field. To generate a 3D random field from a set of Voronoi surfaces, 
we need to thicken the surfaces and make them into fuzzy walls in some way. We take a following procedure in our 
example: first, we randomly distribute points of S in the cubic box. Second, we identify the two nearest random 
points of S from a given point x in the box, and calculate the distances d\ to the nearest random point and d2 to 
the second nearest point. These distances are the functions of x, and they are defined without ambiguity even if the 
number of first and/or second nearest points of S is plural. With those distances, we generate a field p(x) defined by 



where a thickness parameter A is an arbitrary length scale, which corresponds to thickness of the fuzzy walls. In a 
limit of A — > 0, the field value is non-zero only on the Voronoi surfaces. For a finite value of A, the Voronoi surfaces 
are thickened and we have fuzzy walls with Gaussian-like profiles. Below, the fields generated by Eq. (|118|l are called 
Voronoi tessellation fields. 

There are two parameters for Voronoi tessellation fields defined above, i.e., the number density n r of the random 
points initially set to define Voronoi tessellations, and the thickness parameter A. If the thickness parameter is much 
larger than the mean separation of the random points, the generated walls no longer have planar shapes. We fix the 
thickness parameter at one-tenth of the mean separation of points, A = 0.1 n r _1//3 . Thus, we have only one parameter, 
the number density of initial random points. 

In Fig. |21 gray-scale images of 2-dimensional slices of generated 3-dimensional fields are plotted. Four panels have 



A. 



Voronoi tessellation fields 




(118) 
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different number densities of the initial random points. It should be noted that slices of the 3-dimensional Voronoi 
tessellation fields are not 2-dimcnsional Voronoi tessellation fields. We calculate the spatial mean p of the generated 
field and analyze the shifted field / = p — p. 

First, the one-point probability distribution function of a Fourier modulus is calculated. Instead of considering 
many realizations, we consider the distribution of Fourier modulus in a single realization. We take all the modes in 
which absolute lengths of wavevectors are in a certain range, and calculate the distribution of Fourier modulus of 
these modes. Since the field does not have any statistically preferred direction, the function q( 2n \k) of Eq. I|69(l does 
not depend on the direction of k. Therefore, averaging over directions of wavevector does not change the one-point 
distribution function of Fourier modulus. If the function </ 2ra ) (fe) does not change much in a certain range of absolute 
length of wavevectors, the one-point distribution function of Fourier modes in a single realization is still given by 
Eq. I|72|) . In Fig- El the one-point distribution functions of Fourier modulus are plotted. Three cases with different 
numbers _/V r of initial random points of Voronoi tessellations are considered in this figure. The numerical results are 
compared with the analytic formula of Eq. (|66|l . The trispectrum which we need in plotting the analytic curve, 
is numerically calculated from each realization. 

The agreement between the numerical results and analytic formula is quite well. The deviation from the Gaussianity 
is larger for larger N T . When the number N r is small, the scale of characteristic clustering pattern is large, and the 
non-Gaussian correction factor /V of the corresponding scale becomes large. In this example, the lowest-order 
correction works quite well in each case. Even in the case of N r = 10 3 , in which the clustering pattern is quite irregular 
as seen in the upper right panel of Fig. the lowest-order correction is sufficient to describe the non-Gaussianity. 
When the number V r becomes large, deviations from the Gaussian distribution becomes small. This is manifestation 
of the Theorem |21 since increasing the number N T corresponds to increasing the box size relative to the characteristic 
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FIG. 3: The one-point distribution function of normalized Fourier modulus with wavevectors k which absolute magnitudes 
are in a range of |fe| = [15, 25] in units of the fundamental wavenumber 2-k/L. The number of initial random points of Voronoi 
tessellations are N T = 20 3 , 15 3 , and 10 3 from left to right panels, respectively. In top panels, the distribution functions are 
shown. Symbols represent results of numerical realizations. Error bars correspond to Poisson errors of counts. Dotted lines 
represent the Gaussian distribution. Solid lines include the lowest-order correction of non-Gaussianity. The degree of the non- 
Gaussianity /V is indicated in each panel. The middle panels show differences from Gaussian distributions. The bottom 
panels show ratios to the Gaussian distributions. 



scales of the structure. 

Second, we consider the three-point distributions of phases. When there are three wavevectors satisfying fci+fc2 = fc.3 
and there is not any other relations among the three, the three-point phase distribution is given by Eq. i|l(Jl|) . When 
the field is statistically isotropic, the bispectrum should be rotationally invariant. Therefore, averaging over directions 
of wavevector triplets does not change the three-point distribution function of phases. If the normalized bispectrum 
p*- 3 ) does not change much in a certain range of wavevector triplets, the three-point distribution function of phases in 
a single realization is still given by Eq. (|101fl ■ We first consider a wavevector fci , the length of which is in a certain 
range [A;™ m , fc 1 Ilax ]. Next we consider a wavevector k-2, the length of which is in a certain range [£;™ ln , fc™ ax ]. The 
angle 612 between fci and &2 should also be in a certain range [6*™", 0™ 2 ax }- The third wavevector is automatically 
determined by k% = k% + k%. We find all the triplets in a realization according this procedure, and calculate the 
distribution of the phase closure 9k t + 9k 2 — 9k 1 +k 2 - 

In Fig. 01 the distributions of phase closure are plotted. In these examples, the three wavevectors are not aligned, 
and the Eq. (|101|) is applicable. Configurations of the three wavevectors, fci, £2, — &3 = — k\ — k 2 , are approximately 
equilateral triangles. The lowest-order and next-order corrections are indicated by dashed and solid lines, respectively. 
The non-Gaussian correction terms are quite small in each case, and the distribution is almost homogeneous (note 
the scales of vertical axes). 

Similarly, we consider a phase closure of 9k 1 + 9k 2 + $fc 3 — 9k 1 +k 2 +k 3 - Certain ranges for lengths of wavevectors 
fci, fc2, fc3, and their mutual angles are considered, and we calculate distributions of the phase closure of quadruplets 
in these ranges. In Fig. [5] the distributions are plotted. The four wavevectors are not aligned, and the Eq. I|109l) is 
applicable in these examples. The distributions are very accurately homogeneous. 

In the Voronoi Tessellation fields, distribution functions of Fourier modulus can be significantly different from 
Gaussian, while the phase closures are almost homogeneous even though the distributions are strongly non-Gaussian. 
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FIG. 4: Relative distributions of the phase closure 8k 1 + 0k 2 — @k 1 +k 2 - The distributions are calculated for wavevectors 
satisfying |fci| = [19, 21], | fca | = [19, 21] in units of the fundamental wavenumber, and #12 = [118°, 122°], where 812 is an angle 
between the two wavevectors. The number of initial random points of Voronoi tessellations N r , and parameters /VV are 
indicated in each panel. Symbols represent the distributions in numerical realizations. Dotted lines represent the Gaussian 
distribution (i.e., constant). Dashed and solid lines represent the first-order and second-order corrections, respectively, where 
dashed lines are indistinguishable from solid lines. 



B. Lognormal fields 

Our second example of non-Gaussian fields are lognormal fields. A lognormal field is generated from a random 
Gaussian field by an exponential mapping |23j]. First we consider a random Gaussian field (p(x) which has zero mean 
(<f>(x)) — and unit variance ([(f>(x)} 2 ) = 1, where (• • • ) represents a spatial average. This field can have a spatial 
correlation £$(x — x') — (<f>(x)<f>(x')) . From this field, we consider a lognormal field 

p(x) = e 9 ^ x \ (119) 

where g is an arbitrary parameter. Since the distribution of the field 4>(x) is Gaussian, the spatial mean of the 
generated field can be analytically calculated to give (p(x)) = e 9 I" 1 = p. We define a shifted lognormal field 

f{x) = P{X) ~P = «P[g^)-W2]-l (12Q) 
9P 9 

In a limit g — > 0, the above field reduces to a random Gaussian field. In this sense, the parameter g controls the non- 
Gaussianity of the distribution. The non-Gaussianity is large when g is large. We numerically generate a realization 
of random Gaussian field from a given initial power spectrum P^(fc), and a realization of a lognormal field is obtained 
by the above equation. In the examples below, we consider a power-law power spectrum with a high-pass filter 
P<f,(k) oc k n e~ k A / 2 , where we set a spectral index n = and a smoothing scale A = 0.03£ (L is a box size). In 
Fig. El a gray-scale image of 2-dimensional slice of the generated 3-dimensional random Gaussian field <p(x) is shown. 
In Fig. [3 gray-scale images of 2-dimensional slices of the generated 3-dimensional lognormal fields are plotted. The 
lognormal mapping of Eq. (|119fl enhances rare peaks of the initial Gaussian field. For cases the parameter g is of 
order unity or larger, a few peaks are prominent, and other structures are diminished. 
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FIG. 5: Relative distributions of the phase closure 9k 1 +8k 2 +#fc 3 —dk 1 +k 2 +k 3 - The distributions are calculated for wavevectors 
satisfying |fci| = [19,21], \k 2 \ = [19,21], |fe 3 | = [19,21] in units of the fundamental wavenumber, and 6 12 = [98°, 102°], 
023 = [98°, 102°], = [98°, 102°], where 6ij is an angle between the two wavevectors ki and kj. The number of initial random 
points of Voronoi tessellations N T , and parameters /V are indicated in each panel. Symbols and lines are the same as Fig.0] 

The one-point distribution functions of a Fourier modulus are similarly calculated as in the previous case. The 
results are plotted in Fig. [S] The distribution for the case g = 2 is almost Gaussian. We calculate the cases g < 2, 
and they are also indistinguishable from Gaussian. Although the field distribution for the case g = 2 is very different 
from Gaussian as seen in Fig. [3 the distribution of Fourier modulus is almost Gaussian. Therefore, the Theorem[21is 
efficiently achieved in this case. 

In Fig. the distributions of phase closure of three modes, 8^ + 8k 2 ~ dk 1 +k 2 , arc plotted. Choices of three 
wavevectors are the same as in Fig.0] The phase closures exhibit significant deviations from homogeneous distribution 
for g > 1.5. The first-order corrections are sufficient for most of the cases, except when the non-Gaussianity is very 
strong such as g = 3, 4 cases. 

In Fig. I1UI the distributions of phase closure of four modes, 0^ + 0k 2 + &k 3 — 8ki+k 2 +k 3 i are plotted. There are 
again deviations from homogeneous distributions. The degree of the deviations is much less than the phase closures 
of three modes, because the order of non-Gaussian correction terms is second. 

In the lognormal fields, distribution functions of Fourier modulus can be different from Gaussian when the non- 
Gaussianity is strong. Moreover, the phase closures can be inhomogeneous when the non-Gaussianity is strong. 



C. Quadratic fields 



Our last example of non-Gaussian fields are quadratic fields, which are defined below. A lognormal field in the 
previous subsection is generated from a random Gaussian field by an exponential mapping. We similarly define a 
quadratic field, but by a quadratic mapping: 



f(x) = cj ) (x) + h[^(x)-l] , 



(121) 



where <f>(x) is an random Gaussian field which satisfy (4>(x)) — 0, (4> 2 (x)) = 1, and h is an arbitrary parameter. 
We numerically generate a realization of random Gaussian field just in the same way as in the lognormal case, using 
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FIG. 7: Gray-scale images of 2-dimensional slices of the 3-dimensional lognormal fields. Each panel has different non- 
Gaussianity parameter g. Upper left: g = 0.5, upper right: g — 1 .0, lower left: g = 1.5, lower right: g = 2.0. 
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FIG. 8: Same as Fig. [3] but for the lognormal fields. 
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FIG. 9: Same as Fig. 2] but for lognormal fields. It should be noted that scales of the horizontal axes are different. 
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FIG. 10: Same as Fig. but for lognormal fields. It should be noted that scales of the horizontal axes are different. 

the same input power spectrum P<f,(k) oc k n e~ k A / 2 , where n = and A = 0.03L. The parameter h controls the 
non-Gaussianity. For sufficiently small h, the field is essentially random Gaussian. For sufficiently large h, the field 
approaches to a pure quadratic one, / oc 4> 2 — 1. 

In Fig. gray-scale images of 2-dimensional slices of the generated 3-dimensional lognormal fields are plotted. 
The quadratic mapping of Eq. (|121|) enhances peaks of the initial Gaussian field. However, the enhancement is not 
so strong as that in lognormal mappings. 

The one-point distribution functions of a Fourier modulus are calculated, and the results are plotted in Fig. 1121 
The distributions for all the cases are almost Gaussian, in spite of the non-Gaussianity in the fields. Theorem [3 is 
efficiently achieved also in this case. 

In Fig. El the distributions of phase closure of three modes, Bk x + &k 2 ~ #fci+fc 2 > arc plotted. Choices of three 
wavevectors are the same as in Fig. 0] The phase closures exhibit small deviations from homogeneous distribution. 

In Fig. 1141 the distributions of phase closure of four modes, 9k t + 0k 2 +0k a —&k 1 +k 2 +k 3 , are plotted. The distributions 
are almost homogeneous. 

In the quadratic fields, distribution functions of both Fourier modulus and phase closures are similar to random 
Gaussian fields. This result is partly because the volume V is much larger than the characteristic scales of clustering 
in our example. When the volume is effectively small, small deviations seen in Fig ll3l for example, will be enhanced. 

VII. CONCLUSIONS 

In this paper, statistical behavior of Fourier modes in general non-Gaussian fields is studied by explicitly deriving the 
joint distribution function of the Fourier coefficients for the first time. A distribution function for a random Gaussian 
field is very simple: each Fourier coefficients are independently distributed and Gaussian. In a non-Gaussian field, 
there are complex couplings among all modes, and we provide a general framework of entangling this complexity. The 
distribution function is generally expanded by a hierarchy of polyspectra in Eq. I)47p. which has full information of 
statistical properties of the mode couplings. 

The distribution function is formally considered as a series expansion of V -1 / 2 , where V is a total volume of the 
field. If we take a sufficiently large volume, the joint distribution for a particular set of modes approaches to the 
Gaussian distribution. However, this does not mean that the field itself is Gaussian. Non-Gaussianity enters the 
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FIG. 11: 2-dimensional slices of the 3-dimensional quadratic fields. Each panel has different non-Gaussianity parameter h. 
Upper left: h — 0.2, upper right: h = 0.5, lower left: h = 1.0, lower right: h = 2.0. 

distribution function in a volume-dependent way. This is a reason why there appears a total volume in relations to 
define polyspectra by Fourier coefficients. 

The general distribution function is explicitly calculated up to second order of V~ x l 2 , which is given in Appendix lAl 
The distribution function up to this order depends only on bispectra and trispectra. Information of higher-order 
polyspectra is contained in higher-order terms of V~ x l 2 . 

We derive iV-point distribution functions from the general equations. A closed form of expression for the one-point 
distribution function [Eq. (|72|l ] is derived using all the hierarchy of collapsed polyspectra, q^ 2n \k). As a consequence, 
one-point distribution of Fourier phase is shown to always be homogeneous for a random field in a spatially homo- 
geneous space (Theorem As another consequence, the probability distribution function of a particular mode is 
Gaussian in a large- volume limit (Theorem 

For higher-point distribution functions, contributions from lower-point functions are separated, and reduced Ap- 
point distribution functions are introduced. Explicit equations for the reduced functions up to second order are 
given for A^ = 2,3,4. Structures of mode couplings in non-Gaussian fields are generally given for the first time in 
an analytically tractable way. These equations of the general joint distribution function and AT-point distribution 
functions in terms of polyspectra are fundamental equations that can be used to investigate the statistics of Fourier 
modes in non-Gaussian fields, in general. 

The statistics of phase correlations are focused on as an application of the general results. The structure of the 
phase correlations in non-Gaussian fields has been a long-standing issue. We believe our analysis in this work also 
provide a breakthrough in this respect. The distribution function of Fourier phases are straightforwardly obtained 
from our general results. We derive analytic expressions for A-point distribution functions of phases in terms of 
polyspectra. Explicit expressions up to second order are given. 

Regarding the phase correlations, we obtain several theorems which are proven by using the full expression of the 
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FIG. 12: Same as Fig. [3J but for the lognormal fields. It should be noted that scales of the horizontal axes are different. 
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FIG. 13: Same as Fig. [I] but for lognormal fields. It should be noted that scales of the horizontal axes are different. 
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FIG. 14: Same as Fig. but for lognormal fields. 



joint distribution function, i.e., the results without assuming truncated expression of the series. The essence of the 
theorems is that phase correlations among N modes fei , . . . k^ are present only when there is a relation 



- ]Nk N = 0, 



(122) 



where ji, ■ ■ ■ are integers and at least two of them is not zero. Since the zero mode k = is excluded, there 
should be at least one linear relation among the wavevectors with integral coefficients. As a result, there should not be 
any inhomogeneous distribution of a particular phase, because Eq. (|122fl with N = 1 is not possible (This is another 
aspect of Theorem^. For the two-point phase distributions, there are phase correlations only between modes whose 
wavevectors are parallel to each other and corresponding proportional factor should be a rational number (Theorem^J . 

Some results beyond the second order are obtained for phase correlations. A lowest-order contribution to the two- 
point phase distributions is analytically given by Eq. (|1()()|> , even when the proportional factor of the two wavevectors 
are not simple and correlations are of arbitrarily higher order. Similar expressions for V-point distribution functions 
in a limited case that there is only one linear relation of Eq. (|122|l among wavevectors [Eq. I|117|) ]. When there are 
more than one linear relations among the wavevectors, the expression could be more tedious as in Eq. (|112fl . for 
example. 

As numerical checks of the derived equations, we compare some of the analytic equations with numerical realizations 
of three types of non-Gaussian random fields. We consider the Voronoi tessellation fields, the lognormal fields, and 
the quadratic fields. The statistics of Fourier modes differently deviates from Gaussian distributions, depending on 
which type of non-Gaussian field is analyzed. The distributions of Fourier modulus are easily distorted in Voronoi 
tessellation fields, while that in quadratic fields are almost Gaussian. In the three- and four-point phase distributions, 
deviations from Gaussianity in Voronoi tessellation fields and in quadratic fields are small, while that of lognormal 
fields can be large. Although these tendency is not general and depends on our choice of scales and configurations of 
wavevectors, the deviations from Gaussianity appear quite differently from fields to fields. 

The derived analytic results describe the numerical results very well. As for the simple distribution functions 
considered in this work, lowest-order results are sufficient to describe the numerical results in most of the cases. This 
is partly because the volume V is sufficiently larger than the scales of modes we analyze. 

This work provides a basic framework toward understanding the statistical nature of the Fourier analysis. In general, 
the joint distribution function of random variables has the full information on statistical properties of the variables. 
Therefore, all the statistical information on the Fourier modes are contained in the derived joint distribution function. 
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We have analyzed limited number of consequences from this function. Various aspects of the mode couplings in 
non-Gaussian fields, which arises as dynamical effects in physical situations, might be investigated in future work. 
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APPENDIX A: GENERAL JOINT DISTRIBUTION FUNCTION OF FOURIER COEFFICIENTS UP 

TO SECOND ORDER 

In this appendix, the joint probability distribution function of Fourier modes in terms of amplitude A k and phase 
9 k in non-Gaussian fields up to second order is expressed in terms of summations of independent modes. The form of 
Eq. (|57|l is not useful because the modes in each summation with different labels can be the same as described in the 
main text. We expand each summation in Eq. I|57|l so that the different labels refer different modes. The calculation 
requires careful classification of overlapping wavevectors in each summation. After tedious manipulations, the result 
has a form, 
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4A fcl 2 A fc2 A fc3 A fc4 2 cos (0 fel + 6> fc2 - fe J cos (0 fel - fe3 + fe J 



Cl,K 2 ,R 3 ,R4 



- 2^ fcl 2 A fe2 A fe3 cos (20 fel + 9 k2 - fe3 ) - 2A fe2 Afe 3 A fe4 2 cos (0 fe2 + fe3 - 20 fe 

x p (3) (fci,fe 2 , -fc 4 )p (3) (fci, -fes, fe 4 ) , 



uhs 



Y' A kl 2 A k2 A k3 A ki A k5 cos (20 fel - fc4 ) cos (0 fc2 + fe3 - fes ) 

x p (3) (fci,fci,-fc 4 )p (3) (fe 2 ,fe 3 ,-fe 5 ) , 



uhs 

E 

fel ,fc 2 ,fe3,fe4,fe5 



2A fcl A k2 A k3 A ki A k 2 cos (0 fcl - 9 k3 - fes ) cos (0 fc2 - 0fe 4 + 0fe 5 ) 

- A kl A k . 2 A k3 A ki cos (#fci + ^fc 2 - °k 3 - 0fe 4 ) 

x p (3) (fei, -fe 3 , -fe 5 ) p (3) (fe 2 , -fe 4 , fe 5 ) , 



uhs 



Y [ 2A k! A k2 A k3 A ki A k 2 cos (0 fel + fe2 - fc5 ) cos (0 fes - fc4 + fes ) 

fcl,fc 2 ,K 3 ,K4,K 5 

— A kl A k2 A k3 A ki cos (0^1 + 0fc 2 + 0fc 3 — ^fe 4 

x p (3) (fei, fe 2 , -fe 5 )p (3) (fes, -fe 4 , fes) , 



uhs 

E' 

fel ,fe 2 ,fe 3 ,fe4,fe5 



-A fel A fc2 A fe3J 4fc 4 Afc 5 2 cos (0 fcl + 0fc 2 - 6<fe 5 ) cos (0 K3 + 6» fc4 - fes ) 



-A kl A k2 A k3 A k4 cos (0fcj + 0fe 2 — 0fc 3 — ki ) 



p (3) (fci,fc2,-fe 5 )p (3) (fe 3 ,fe 4 ,-fe5), 



uhs 



. V 1 

-A kl A k2 A k . 3 A ki A k5 A k6 cos (0 fel + 9 k2 - # K3 ) cos (6 ki + 9 ks - 6 k(i ) 

x p (3) (fei, fe 2 , -fc 3 )p (3) (fe 4 , fes, -fee) , 



fcl ,k2 ,&3 ,^4 ,&5 ,^6 
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In the above equations, the abbreviated symbol s k indicates that any two of wavevectors fci, k2, ■ ■ ■ appeared 

(i) 

in the summation should be mutually different. It is important to note that each quantity vanishes when any 
one of the modes appearing in a summation is integrated over, because of Eqs. i|til[l - Ht)3[) . 
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